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Abstract 



> 

O 

^^ I The dynamics of spin-boson systems at very low temperatures has been 

^^ I studied using a real-time path-integral simulation technique which combines 

^ ' a stochastic Monte Carlo sampling over the quantum fluctuations with an 

^' 

O '_ exact treatment of the quasiclassical degrees of freedoms. To a large degree, 

this special technique circumvents the dynamical sign problem and allows 

0>\ I the dynamics to be studied directly up to long real times in a numerically 

o ' exact manner. This method has been applied to two important problems: 



(1) crossover from nonadiabatic to adiabatic behavior in electron transfer 
reactions, (2) the zero-temperature dynamics in the antiferromagnetic Kondo 
region 1/2 < K < 1 where K is Kondo's parameter. 
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I. INTRODUCTION 

Spin-boson systems provide archetypical models for many low-temperature dissipative 
quantum tunneling systems |[l|-^. Applications include diverse problems such as the ob- 
servability of macroscopic quantum coherence in a superconducting quantum interference 
device (SQUID) |^, interstitial tunneling of light particles such as hydrogen in metals 0, 
and many others discussed in Ref. |T[]. Two other examples are of special importance to this 
work. First, the so-called Kondo problem concerns localized spin impurities in nonmagnetic 
metals 0. The second comes from the realm of chemical physics — in certain parameter 
regions, the spin-boson model provides a generalization of the Marcus model for electron 
transfer reactions 0-^1 . The wide range of applicability of the spin-boson model stems 
from its apparent generality - by coupling two exactly solvable models bilinearly, namely a 
two-state system (spin) and an infinite-dimensional harmonic oscillator (boson) bath, one 
obtains a nontrivial description for many dissipative systems [|10 . 



The dynamics of the spin-boson system has been widely studied using the Feynman- 
Vernon influence functional method [|ll|,|l^, mostly in conjunction with instanton techniques 



0. Much of the physics of this model has been unraveled by analytical methods, though 
an exact solution is not possible (except for some special parameter values). One particu- 
larly useful analytical approximation is the noninteracting blip approximation (NIBA) [Q. 
Although the NIBA has been rather successful, there are important regions in parameter 
space where this approximation is expected to fail. 

Very recently, numerical techniques for computing the dynamics of dissipative two-state 
systems have been developed [p!3| -[T6|. These real-time quantum Monte Carlo (QMC) simu- 
lations have confirmed the NIBA predictions quantitatively in many cases, and in addition 
have provided information in regions where the NIBA breaks down. However, due to the 
fundamental dynamical sign problem inherent in these calculations [p!7|-p2|, the numerical 
computations are restricted in their accessible range of real times. The dynamical sign prob- 
lem arises because at long times a large number of interfering paths contribute, leading to a 



very small signal-to-noise ratio. In effect, the simulation becomes unstable. Consequently, 
even numerically exact QMC methods have not been able to resolve many important ques- 
tions concerning the behavior of the spin-boson system, especially at low temperatures and 
long times. 

It is important to point out that many other approaches to dealing with the dynamical 
sign problem have been developed during the last decade. For example, methods based on 
related filtering techniques |[17| , [T8[| , optimized reference systems fl^, or analytic continua- 



tion procedures following a conventional imaginary-time QMC simulation pO[ have been 
proposed. Similar to the one described in this work, the first two methods attempt a direct 
simulation in real time. While we restrict our attention to discrete (tight-binding) systems 
which is adequate for many problems in low-temperature physics, these previous methods 
are designed for extended systems. On the other hand, the third method makes use of an 
imaginary-time simulation that has no dynamical sign problem to obtain numerical data 
which are then analytically continued to real times. To date this procedure has only be 
used to find (e.g., electronic) spectral densities. While real-time information follows from 
the knowledge of these densities, such calculations have yet to be performed and hence the 
effectiveness of this approach to real-time dynamics is untested. We also point out that a 



similar numerical problem exists for QMC simulations of fermionic many-body systems p3 
This fermion sign problem has a different origin from the one dealt with here. It is due to 
the antisymmetry of fermionic wavefunctions, and makes even imaginary-time simulations 
problematic. If simulated in imaginary time, however, the spin-boson model does not pose 
any sign problem. 

In this paper we propose a new simulation method to study the spin-boson dynamics 
numerically. Like our earlier approaches, this technique is based on a discretized path 
integral representation of the dynamical quantities of interest, but differs from them in 
that it exploits the symmetry properties of the influence functional. This allows for an 
exact treatment of the numerically problematic quasiclassical paths, and one is left with 
a stochastic Monte Carlo sampling of the quantum fluctuations alone. Related methods 



have been used previously to study quantum Brownian motion pT|j2^ , and we have used a 
similar algorithm to study the primary electron transfer steps in the bacterial photosynthetic 
reaction center |^. Furthermore, this idea of exploiting the symmetries of the influence 
functional has also led to an efficient simulation method for computing the mobility and 
diffusion coefficient of a dissipative particle in an infinite tight-binding lattice [^. The 
algorithm for this diffusion problem, however, is rather different from the one employed in 
the case of a system with a small number of tight-binding states. In this article, we provide 
a detailed description of our simulation method for the dissipative two-state system; the 
necessary generalizations for the case of more than two states should be obvious. 

We apply this special method to the spin-boson model in the low-temperature and the 
small bath cutoff regions. The accuracy of the NIBA in some other parameter regions has 
been confirmed by our earlier simulations already, and the numerical results presented in 
this work focus on these previously inaccessible regions. Since the present algorithm is 
more powerful than previous methods, we are able to study longer times and much lower 
temperatures. 

Data are presented for two practical problems of current interest: (1) Electron transfer 
in a condensed phase environment, where we have examined the transition from nonadia- 
batic to adiabatic behavior for both the high-temperature and the low-temperature case. 
The high-temperature results should coincide with classical Marcus theory, whereas for low 
temperatures we would expect quantum effects commonly attributed to nuclear tunneling 
0. (2) The dissipative two-state system is also significant for the antiferromagnetic Kondo 
problem of a localized impurity embedded in a nonmagnetic metal [Q. The anisotropic 
Kondo Hamiltonian is related to the spin-boson model in a certain parameter region, and 
we have studied the most interesting antiferromagnetic case 1/2 < K < 1, where K is 
Kondo's parameter, at zero temperature. An important question to be addressed in this 
parameter region relates to the destruction of quantum coherence. For K < 1/2, one can 
observe quantum coherence (oscillatory behavior) in the zero-temperature dynamics of the 
dissipative system. The NIBA predicts such coherent behavior to be completely destroyed 



for K > 1/2. Unfortunately, the justification for NIBA is suspect in this special parameter 
region, and our method offers an unique way of studying this region. 

In Section II we present our simulation method for the dissipative two-state system. The 
results for the crossover between nonadiabatic and adiabatic electron transfer are discussed 
in Section III, followed by a discussion of the dynamics in the antiferromagnetic Kondo 
region at very low temperatures. Some final remarks and conclusions are given in Section 
IV. 

II. SIMULATION METHOD 

In the following we present a dynamical simulation technique for dissipative tight-binding 
models with a finite number of states. For simplicity, the discussion is restricted to the case 
of two states, which leads to the often-studied spin-boson Hamiltonian p| 

H = Ho + Hb + Hi (2.1) 

Pa I 1,^ ,2 / ^ ^a 



-{nA/2)a., + {he/2)a, + J2 






2ma \ maUJ^ 

The parameters in the free Hamiltonian Ho describing the isolated two-state system are 
the tunnel matrix element A — in the parlance of electron transfer theory, A/2 is the 
electronic coupling between the two different redox sites, the external bias e corresponds to 
an asymmetry between the two localized energy levels, and a^ and a^ are the usual Pauli 
matrices. The bath is described by harmonic oscillators {xa} which are bilinearly coupled 
to the spin operator o"^. This type of coupling is reasonable for the problems considered in 
this work (and for many others); for a justification, see, e.g., Refs. [|l]-3]. 

Within this model, the bath parameters enter only via a single function called the spectral 
density 

J{u;) = ^Y.-^^i^-^c.), (2.2) 

which should have a continuous form in the limit of infinitely many bath oscillators. The 
spectral density then determines the bath correlation function |^ 



nh Jo smh.[ujhp/2] 

which is defined for complex-valued times z = t — ir. The distance between the localized 
states is given by a lengthscale a, and (3 = l/ksT. Here, we limit our attention to the case 
of an Ohmic spectral density which has the form 

J{uj) = {27ThK/a^) u e"'"/^^ . (2.4) 

This spectral density has a characteristic low-frequency behavior J{uj) ~ rjuj, where rj is 
the usual Ohmic viscosity. The system-bath coupling strength is measured in terms of the 
dimensionless Kondo parameter K. The timescale distribution of bath motions is described 
by a cutoff frequency Uc- For many problems in low-temperature physics, this cutoff fre- 
quency is taken to be the largest frequency scale in the problem. In the case of electron 
transfer, the same spectral density with some intermediate value for Uc is most appropriate 



for a realistic description of many polar solvents p6 . 

Two dynamical quantities of interest for this model are the symmetrized time correlation 
function 

C(t)=Re(a,(0)a,(t)) 

= Z-iReTr (e-^^a^e^^^/V^e"*^*/") (2.5) 

with Z = Tt e~^^ , and the time-dependent occupation probabilities for the two sites P+{t) 
and P- (t) which can be expressed in terms of a single function 

P{t) = P+(t) - P4t) = (e'^*/V,e-'^*/^) (2.6) 



with the initial condition P(0) = 1. From a comparison of Eqs.(2^) and (p.6|) we observe 



that the two quantities differ only by the way the system is prepared initially. For P{t), the 
system is held fixed in the + state until t = with the bath being unobserved (factorized 
initial condition). On the other hand, the more realistic situation for many experiments is 
represented by equilibrium states of the total system at t = as described by C(t). 
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To numerically compute the two-state dynamics, we employ a discretized path-integral 
representation of the dynamical quantities |]l2|-p!5|]. The correlation function {crz{^)<^z{t)) can 
be regarded as the probability amplitude for a sequence of steps in the complex-time plane. 
In particular, one propagates along the Kadanoff-Baym contour 7: 2 = 0^t^0^ —ihp, 
and measures cr^ at z = and z = t. Of course, there are several other possible choices of 
this contour due to the cyclic structure of the trace. 

To parametrize the paths, we use q uniformly spaced discrete points for each of the 
two real-time paths and r points for the imaginary-time path (see Fig.|lD. Hence, there are 
N = 2q + r points in total. The time discretizations are 



A. 



t/q , l<j<q 

-t/q , q+l<3 <2q (2.7) 

-ihf3/r , 2q + l<j<2q + r = N, 



and the complex time after i steps is Zi = Y,]=i ^j- The construction of the discretized path 
integral then proceeds by inserting complete sets 

= Jdr,\r,){ri\ (2.8) 

at each discretization point Zi {i = 2, . . . , N). The vector Tj represents the state ctj = ±1 of 
the two- level system as well as the environmental degrees of freedom {xq^j}. 

To disentangle the short-time propagator, we use a (symmetrized) Trotter formula, 
exp{-iHAj/n) = exp{-tH'Aj/2h) expi-tHoAj/h) exp{-iH'Aj/2h) + 0{A][Hq, H']) , 

(2.9) 



where H' = H — Hq is diagonal in the representation (p.8|) . The free part Hq of the 
Hamiltonian ( p.l|) leads to the short-time propagator 



K{aj,aj+i) = {(rj+i\exp{-iAjHo/h)\aj) (2.10) 

which can be evaluated exactly. Thereby we arrive at 



M^)--^{t')) = ^^,^^^,t,,...,^, , (2.11) 

which converges to the true path integral as the number of discretizations N ^ oo. The 
discretized action S'[ri, . . . , r^v] is a complex- valued sum of the actions picked up in separate 
parts of the contour 7. To compute the correlation function in Eq. (|2.11 ), one can average 



over all pairs of spins {a', a"} separated by a time f along the contour 7. This allows 
us to compute the dynamical quantities for all times tk = kt/q (where k = 0, . . . ,q) from 
one single Metropolis trajectory. Because of the cyclic structure of the trace in Eq.( |2.5| ), 
tn+1 = ri. 

Since the bath is made up of harmonic oscillators, one can integrate out the environmental 
degrees of freedom analytically. After performing this integration, the bath-plus-coupling 
part H' of the Hamiltonian leads to an influence functional $[cr] in terms of the spins {ctj} 
alone [|ll] , |T2| . As a result, the correlation function takes the form 



{crM^zit')) = ^E^xp (-^cr]+J2\nK{a„a,+,)] a'a" , (2.12) 

^ {a} V i I 

where Z = I]{o-} exp(. . .) (the exponent will be referred to as "the action" henceforth). In 
the continuum limit N —^ 00, the nonlocal influence functional is given in terms of the bath 
correlation function L{z) introduced in Eq.( |2.3| ), and one finds with cXi -^ a{z) 



^a] = f dz f dz' a{z)L{z - z')a{z')/4: . (2.13) 

The integrations in the complex-time plane are ordered along the Kadanoff-Baym contour 
7. Note that L{z) fulfills the important symmetry relation 

L{z-ih(3) = L{-z) , (2.14) 

which implies certain symmetry properties of the influence functional. 



In discretized form, the influence functional is given by [12 



N 

'^W] = kJ2 ^jLjkCTk/^, (2.15) 

j,k=i 



with the complex- valued influence matrix Ljk. We observe from Eqs.( p.l^ ) and (|2.15| ) that 



the (discretized) spin-boson problem is isomorphic to a classical one-dimensional Ising model 
with long-range complex- valued interactions [|12],|13|. The influence matrix is the average 
value of the influence functional interaction between two points Zj and Zk 

Ljk = Lkj= dz'j dz'j, L{z'. - z'f.) {ioi j > k) , (2.16) 

where Cj is one discretization on the contour centered at the point Zj, i.e., 

{z G Ci\zi - A,_i/2 <z <z, + Ai/2} . 

The remaining time integrations in Eq.( ^.16D can be carried out easily, and one flnds with 

^jk ^j ^k 

L,k = g(A,fe + (A, + Afc_i)/2) + Q{A,k + (-A,_i - Afe)/2) 

- Q{A,k + (-A,_i + Afc_i)/2) - Q{Ajk + (A, - Afe)/2) , (2.17) 

where Q{z) is the twice- integrated bath correlation function with Q{0) = 0, i.e., 
(fQ{z)/dz'^ = L{z). Of course, this function exhibits the same symmetry property ( p.l4[ ). 
Finally, the diagonal elements are given by 

L,, = 2Q((A,_i + A,)/2) . 



A detailed derivation of the influence matrix can be found in Refs. ||T5|j24 



Since the action for each spin conflguration ("a path") is complex- valued, a stochastic 
Monte Carlo evaluation of the resulting isomorphic Ising chain is faced with the dynamical 
sign problem. In the past, we have partially circumvented this problem either by transform- 
ing to a continuous spin representation and applying stationary-phase Monte Carlo (SPMC) 
methods |]T3| , p3| , or by introducing a local flltering function in discrete state space |1^,|T6[. 
The latter approach is related to ideas like the stationary-phase approximation and allowed 
for a study of many phenomena on an intermediate timescale. 

Here, we observe that a much more efficient method can be constructed when one takes 
into account the symmetry properties of the influence matrix due to Eq.( p.l4| ). These 
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symmetry relations can be expressed mathematically as a set of index relations [we use 
i,j,k = l,...,q + l, and n, m = 2, . . . , r] 

L2q+2-i,2q+2-j = -^jj 

L2g+2-i,j = -L*j for i > i 
L2q+2-i,j = -Lij for i > j 

L2q+2-i,i = — Rc La 
L2q+m,2q+2-j = —L2q+m,j for J > 1 

Li,q+1 = L2q+2-i,q+l = -^2g+m,g+l = , 



which can be proved easily from Eq.( p.l7 ) 



The benefit of exploiting these relations is realized upon switching to a new spin repre- 
sentation. To that purpose, we introduce the sum and difference coordinates of the forward 
((Tj) and backward (cr') real-time spin paths and rename the imaginary-time spins (o"m), 

e. = (^. - ^;)/2 (2.18) 

where a' = (T2q+2-j in the old notation (j = 1, . . . ,q + 1). Thus we first relabel the spins on 
the three pieces of the Kadanoff-Baym contour and then form the said linear combinations. 
A physical meaning can be assigned to these new spins: {rjj} describe the propagation along 
the diagonal of the reduced density matrix and can thus be identified with quasidassical 
paths, while {^j} describe the ofF-diagonaliticity of the reduced density matrix and can be 
identified with quantum fluctuations |2^. According to the definitions ( p.l8| ), the new spins 
can take on three possible values ^,1] = —1,0,1, but they are not entirely independent 
because either C,j or rjj has to be for the same j. 

Written in terms of the new spins and exploiting the index relations above, the influence 
functional takes the form 
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r q 

m,n=2 j,k=l 

q q r 

jyk=l j=l m=2 



+ ^71 I 2^ ^m [L2q+m,l + L 



2q+m,2q+l\ 



\m=2 



The elements of the matrices appearing in Eq. ( |2.19|) are 

^ mn ^2q+m,2q+nl ^ 
^jk = Re Ljk 

Xjk = ImLjfc 

^jm J-'j,2q+ml ^ ■ 

It is worth mentioning that these matrices are real-valued [with the exception of Zjm]- The 
meaning of the five terms in Eq.( |2.19|) is as follows. The first term describes a self- interaction 
within the imaginary-time segment. Due to the second term, the dissipative bath will damp 
out the quantum fluctuations, and the system is likely to be found in a diagonal state 
characterized by ^ = 0. This point will later be important with regard to the choice of a 
suitable Monte Carlo weight. The third term is a bilinear interaction between quasi classical 
paths and quantum fluctuations, and the fourth term describes a similar interaction between 
the imaginary-time spins and the quantum fluctuations. The last term is a preparation 
term, coupling rji to the imaginary-time spins. Remarkably, there is no self-interaction in 
the quasiclassical paths, and they are coupled only linearly to other degrees of freedom. 

This observation is crucial for our computational procedure, since it allows for an exact 
treatment of the quasiclassical paths. For any given quantum fluctuation path, the path 
summation over all allowed quasiclassical paths can he carried out in an exact manner. To 
elucidate this, we first examine the free action due to Hq. The imaginary-time contribution 
can be put into the matrix elements Ymn by simply adding | lntanh(fi/9A/2r) to Ym,m+i and 
Ym+i,m', in case an external bias e is present, the action has to include an additional term 



{h(3e/2r) J2m '^■m HTH • Regarding the real-time paths, we proceed in a different way. For the 
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isolated two-state system, we have (in terms of the original spins) a product of the form 
[cf. Eq.(^)] 

ni^(a,,a,+i)i^*(^>;+i)- (2.20) 

i=i 

If we now switch to the {C,,f]} representation and perform the summation over all t] spins 

(while keeping the ^ configuration frozen), we obtain a matrix product for Eq. (|2.2CI|) . Of 

course, one has to account for the ?7-dependent terms in the influence functional (p.l9|) 



during this procedure. In the end, the complex-valued contribution of all these terms for a 
given C, configuration takes the form 

J[^]= E E (m|V«...V('')|r^,+i), (2.21) 

r?i=0,±l»7,+i=0,±l 



where the (3 x 3) matrices V*^''^[^] are defined by |^ 



g 
k>j 



(2.22) 



Each of the matrices V*--^^ depends on all ^k spins with k > j; however, the "free" part 
K X K* is determined by ^j and ^j+i alone. Clearly, J'[^] can be evaluated with a simple 
matrix multiplication routine, leading to a numerically exact and efficient treatment of 
the quasiclassical paths. Note that the remaining part of the influence functional is real- 
valued — with the exception of the fourth term in Eq.( |2.19| ), which is generally very small 



— indicating that much of the dynamical sign problem has been relieved by treating the 
numerically problematic quasiclassical paths in an exact manner. 

In effect, the factor J'[^] contains all contributions from Hq and, in addition, the third 
and fifth term of the infiuence functional ( |2.19| ). The correlation function can thus be written 



as 



1 I '" - - 

{5}=-l,0,l{5-}=±l V m,n=2 

q q r \ 



2 E ^j^jkik ^YiZ^ ijZjmCTrr,. I CT' (T" 



j,k=l j=l m=2 
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Therefore we are left with the task of summing over the imaginary-time spins {am} and 
the quantum fluctuations {C,j}, which is conveniently done via Monte Carlo (MC) sampling. 
The suitable MC weight for the imaginary-time spins is straightforward, 

^imag[0-] ~ exp -i ^ O-mYmnO-n " Re ^^^Zj^O"^ 

Since the influence functional forces the quantum fluctuations to stay near the diagonal of 
the density matrix, we first try to use the following MC weight for the quantum fluctuations 

The problem with this weight, however, arises for small system-bath couplings where the 
damping of the quantum fluctuations becomes weak. In this case, the importance sampling 
would become very inefficient for small coupling K. Hence, in a next step, we try the product 

where J[E] has been defined in Eq.( p.21] ). This weight function considers both the damping 
of the quantum fiuctuations due to the infiuence functional and the integrated-out quasi- 
classical paths. 

Unfortunately, there is another problem with this weight, similar to the case of the mul- 
tistate algorithm [^,0. This problem arises since for correlation functions like (cr2(0)o"2(t)) 
one has to compute the ratio of two quantities. It turns out that certain spin configurations 
only contribute to the denominator but not to the numerator (and vice versa). Due to 
this exclusivity problem, one will not be able to access all relevant spin configurations {^} 
contributing to the numerator when using W[^] alone as the Monte Carlo weight. We can 
circumvent this problem by observing that for the numerator, one has to compute terms like 
[where ai, 0^2 = ±] 

Ji%M= E E (2-23) 

77i=0,±lr7g+i=0,±l 
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The projection operators -ff+ = (1 + crz)/2 and H 
have the //-representation (for a given S,) 



[1 — (y.^j2 onto the two spin values 



^+(e = o) 



' 1 0^ 








ii±{i = ±1) 






0^ 




1 




1/2 




Finally, an appropriate positive definite Monte Carlo weight function can be constructed 



w^[e]~neai[e] E E pir, 



(fc) 



(2.24) 



Using W [^] as the weight allows us to carry out an efficient Monte Carlo sampling of the ^ 
spins. 

Our QMC algorithm employs single particle Metropolis moves as well as moves which 
allow kinks to translate along the spin chain |]16|. Single-particle moves attempt to change 
one spin C,k = ~1, 0, 1 to a new value, whereas kink moves attempt to change two adjacent 
spins simultaneously. The imaginary-time spins are sampled from Pimag[o"] using single 
particle moves, i.e., one tries to fiip a single spin a^ = ±1- During one MC pass, the 
single-particle moves are attempted once for every spin, and the kink moves are attempted 
for every pair of spins with ^^ 7^ ^fc+i- Typical acceptance ratios for these types of moves are 
~ 15%, and we take samples separated by 5 MC passes. This ensures that the MC samples 
are sufficiently uncorrelated, since roughly half of the spins have been assigned new values 
between two subsequent samples. Numerical results are then obtained from several 10,000 
samples, with statistical errors always below 5% for the data reported here. 

The calculations were carried out on an IBM RISC 6000/580 workstation, at an average 
speed of 2 CPU hours per 10,000 samples (for q ^ 80, r = 0). As mentioned earlier, the 
dynamical quantities of interest can be sampled from a single Monte Carlo trajectory for all 
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times t' <t since the remaining part of the Kadanoff-Baym contour can be integrated out. 
Furthermore, the quantity P{t) can be calculated using the same code by simply removing 
all spins from the imaginary-time path, i.e., by putting r = 0. Finally, to ensure that 
the Trotter error is sufficiently small, we have to keep the discretization numbers q and r 
large enough. This is checked by systematically increasing these numbers until convergence 
is reached. For all results presented here, the Trotter error is negligible compared to the 
statistical errors due to the stochastic MC sampling (which are less than 5%). 

We close this section with some remarks concerning the relation between this method 
and our earlier techniques. The algorithm presented here can be thought of as an optimized 
(but essentially nonlocal) filtering method for discrete-state systems as proposed by Mak 
H^ . The optimization is achieved by exploiting the underlying symmetries of the infiuence 
functional. This makes this method superior in the sense that we can study longer times 
with less computational effort. Since the dynamical sign problem grows exponentially with 
increasing time, most of the results discussed in the next section cannot be obtained by any 
former technique. An approximate measure for the performance of the different real-time 
QMC algorithms for the spin-boson problem can be obtained from the maximum real time 
^max defined as the upper time limit of the respective method. For times t > tmax, the 
large statistical errors caused by the dynamical sign problem (more than k, 20%) will render 
the simulation results useless. This is quantified in Table I. The comparison with earlier 
methods demonstrates that the spin-boson dynamics can now be studied up to much longer 
timescales despite the dynamical sign problem. The gain is most significant at very low 
temperatures, but is also important at higher temperatures. 

III. SIMULATION RESULTS 

In this section, we present dynamical simulation results for two different problems, 
namely the crossover between nonadiabatic and adiabatic electron transfer and the dynam- 
ics in the low-temperature antiferromagnetic Kondo region. We will restrict our attention 
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to the symmetric case e = here, and we also consider only Ohmic spectral densities of the 
form (p.4|). For numerical results in other parameter regions of the spin-boson model, we 



refer to our previous work [0-|T6| 



A. Crossover from nonadiabatic to adiabatic electron transfer 

The spin-boson system is an adequate model for many electron transfer reactions in con- 
densed phase systems 0. The electron transfer rate is in general determined not only by the 
overlap of the electronic wavefunctions localized on the redox states (which is proportional 
to the tunnel splitting A), but also by the properties of the solvent environment. For a 
charge transfer to occur, a specific large-scale reorganization of the solvent is required to 
achieve the resonance condition necessary for electronic tunneling. Using linear response 
theory for a description of the solvent modes and a tight-binding model for the redox states, 
one arrives at the spin-boson model. In our study, we have taken an Ohmic spectral density 
for the bath with a finite cutoff frequency Uc which can be regarded as a free parameter 



The solvent is described by a continuous spectral density peaked around a characteristic 
bath mode frequency tUc, and the classical reorganization energy corresponding to this bath 
is h\ = 2KhuJc, where K is Kondo's parameter. For typical electron transfer reactions, 
the reorganization energy fulfills the condition A/ A ^ 1. The nonadiabatic regime of 
electron transfer is defined by small electronic couplings, A/cJc <^ 1, whereas the adiabatic 
limit corresponds to electronic couplings of the order of Uc (or larger). We note that the 
nonadiabatic limit is realized in most biological and in many chemical systems; nevertheless, 
the adiabatic limit is also important for a description of many chemical electron transfer 
reactions. 

The high-temperature limit of electron transfer is well understood within the framework 
of classical Marcus theory [^ . The rate can be factored into an equilibrium Boltzmann factor 
containing the activation free energy for the required global bath fluctuation, and a Landau- 
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Zener factor for the transition probability once this Landau-Zener region has been reached 
PJ. For symmetric electron transfer reactions, Marcus found a bath activation energy hX/4, 
and combining this with a conventional estimate for the Landau-Zener factor M, one obtains 



a formula for the total (forward plus backward) rate |^ 

which is valid for hf3ujc <^ 1. The Landau-Zener factor contains a frequency scale tUb ~ A 
reminiscent of transition state theory which is usually applied to the adiabatic limit [^ . Note 
that for small electronic couplings, the rate ~ A^ which is just the golden rule behavior. 
In the adiabatic limit of large A, however, the rate becomes independent of A, and the 
dynamics is totally solvent-controlled. We note that for the spectral density (|2.4| ), one can 
obtain analytical expressions in the nonadiabatic golden rule limit for both the high- and 
low-temperature rate (T' = huc/k-B provides a rough measure for the temperature below 
which quantum effects due to nuclear tunneling become important). For the special values 
K = 1/2 and K = 1, the crossover behavior from high to low temperatures can be solved 
exphcitly PB| . 

In Fig.|^ we show some results for P{t) in the high-temperature limit; in this parameter 
region, the effects of the initial preparation are negligible, so C(t) = P(t). For a study of the 
crossover between nonadiabatic and adiabatic behaviors, all model parameters except A are 
kept fixed. In the case of small A, the simulation results exhibit a monoexponential decay 
on the golden rule timescale. However, for larger A, the dynamics becomes progressively 
more complex. After a fast initial transient, the decay slows down; fitting this slower decay 
to an exponential law, one can again extract a rate for this adiabatic situation. 

The A-dependent high-temperature rates measured in units of uJc are plotted in Fig.|^. 
Clearly, the nonadiabatic limit is nicely reproduced by the simulations for small A, and the 
rates ~ A^. In the adiabatic limit, the rate constant is approximately independent of the 
magnitude of the electronic coupling. These results are in agreement with the conventional 
Landau-Zener prediction ( p.l|) in the limit of high temperatures. Remarkably, the data in 
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Fig.^ show that the golden rule formula is accurate for values of A/cUc as large as ~ 1.5. 
Note also that we find a monotonic dependence of the rate on A, in contrast to the findings 



of Skourtis et al. |^ which are based on a rather simplistic Hamiltonian and would predict 
a maximum in the rate as a function of A. 

Next we study the low-temperature region {h(3ujc = 2.5), where the classical rate formula 
( p. ID is not expected to hold. Again, we find complete agreement with the nonadiabatic 
golden rule formalism for small values of A/cUc, with P{t) decaying monoexponentially on 
the golden rule timescale (with a fast initial transient). When increasing A/cJc beyond ~ 0.5, 
however, the dynamics exhibits an oscillatory behavior as shown in Fig.^. This prevents a 
meaningful estimate for the decay constant (but when increasing K to significantly larger 
values, we expect that the dynamics will become totally incoherent again, at least for not 
too small cUc/A). We note that in the limit uj^/A -^ (but A finite), the dynamics can be 
solved exactly. This is the case of a strictly adiabatic bath |^ , where the dynamics is always 
oscillatory. From the simulations we also observe that the nonadiabatic regime is confined 
to increasingly smaller values for A/uc when the temperature is lowered. 

In conclusion, the simulations confirm the classical picture for the crossover from nonadi- 
abatic to adiabatic electron transfer in the high-temperature limit and, in addition, provide 
a dynamical explanation. In the low-temperature limit, however, the dynamics becomes 
oscillatory unless the system-bath coupling is made very large. 

B. Antiferromagnetic Kondo region: Low-temperature dynamics 

We next turn to a different problem. The determination of the dynamical quantities of 
interest in the low-temperature Kondo region characterized by a Kondo parameter 1/2 < 
K < 1 has been a long-standing problem. 

We start with a brief discussion of the relationship of the spin-boson model to the Kondo 
problem 0. The Kondo Hamiltonian in its simplest form describes a spin-^ impurity inter- 
acting with a band of free electrons via isotropic exchange scattering. A particularly useful 



method for studying the low-temperature spin relaxation dynamics of the Kondo problem 
(the equilibrium properties are well understood 0) employs a bosonization procedure to 
map it onto a spin-boson problem with Ohmic dissipation ||^. The case 1/2 < K < 1 then 
corresponds to the interesting antiferromagnetic Kondo problem. The important dynamical 
quantity in the Kondo problem is the imaginary part of the frequency-dependent spin sus- 
ceptibility, ^("{uj). It can be expressed in terms of the Fourier transform of the correlation 
function Q 

C{u) = ncoth{(3huj/2)x"{uj) . (3.2) 

Therefore a computation of C{t) will yield all relevant dynamical quantities (structure factor, 
dynamical susceptibility, etc.) of the Kondo problem. 

Unfortunately, the noninteracting blip approximation (NIBA) by Leggett et al. cannot 
be justified in this region for temperatures below the Kondo temperature defined as Tk = 
hAr/kB, where A, = A{A/uj^)^/^^-'^\ In fact, if one equates C{t) = P{t) (which NIBA 
predicts is true), the NIBA would imply certain unphysical properties of the related Kondo 
problem, such as a divergence in the susceptibility as T -^ 0. 

Furthermore, this parameter region is also interesting in the context of macroscopic 
quantum coherence, since there could be remnants of coherent (oscillatory) behavior in the 
zero-temperature dynamics for 1/2 < K < 1. For this region, NIBA predicts a complete 
destruction of quantum coherence. We have studied P(t) at zero temperature in order to 
check this prediction. 



Some numerical results for this parameter region have already been given in Refs. |Tj,|T5 
and a biexponential behavior has been found. However, due to the limitations of our earlier 
techniques, these simulations were restricted to very small values for the cutoff {uJc/A = 1.25 
has been used in Ref. |]15[). Furthermore, our previous algorithm did not allow for a study of 
extremely low temperatures, and we have considered temperatures only slightly below Tk. 
These shortcomings can be overcome by the algorithm discussed in Sec. II, and we are now 
able to reach both the scaling region tUc/A ^ 1 and the zero-temperature limit for P{t). At 
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zero temperature, the timescale of the dynamics should be set solely by the frequency |^ 



Aeff = [cos{7cK) r(l - 2K)]i/2(i-^) A, (3.3) 

where T{x) denotes the Gamma function. This effective frequency scale is equal to A at i^ = 
0, becomes equal to ttA'^/2uJc for K = 1/2, and shrinks to zero as i^ — ;> 1. Since the cutoff 
Uc should enter the dynamics only via a renormalization of this effective frequency scale, we 
use the dimensionless time y = Acgt; at zero temperature, the only system parameter left is 
the Kondo parameter. As shown by Grabert and Weiss [Q, the NIBA solution for T = 



and K < 1 can be written in terms of the Mittag-LefHer function [^ , 



PNIBA(y)=i?2(l-i.)(V^'-^^)- (3.4) 

In the parameter region 1/2 < K < 1, this describes a purely incoherent relaxation. 

Before discussing the zero-temperature limit for P{y) and 1/2 < K < 1, we first present 
data for the correlation function C{t) at K = 1/2. This special case has been solved 
exactly by Sassetti and Weiss |^, and one finds for sufficiently low temperatures that C{t) 
approaches zero from below at long times. This is reproduced by our simulations shown 
in Fig.Q the cutoff chosen here {uc = 6A) is already large enough to ensure the validity 
of C(Jc/A ^ 1. Note that the corresponding exact solution for a factorizing initial state, 
P(t) = exp(— Acfft), does not exhibit this behavior. 

Finally, in Fig.^ we show QMC results for the exact T = dynamics of P{y) in the 
Kondo region 1/2 < i^ < 1. The cutoff chosen here is within the scaling region cUc/A ^ 1 
so that the dynamics should depend only on the effective timescale y = Aefjt. Indeed, as long 
as uJc/A > 5, we found that a change in the cutoff enters the dynamics solely via Eq.( |3.3| ). 
Since the frequency Aefr becomes extremely small with increasing K, it is not possible to 
numerically study the T = dynamics on timescales of the order A^Jg^ for K larger than 
~ 0.75. Here, we have restricted ourselves to K = 0.6 and K = 0.7, see Fig.^. For K = 0.5, 
the QMC data for P{t) coincide with the exact solution. 

It is not possible to fit the numerical data in Fig.|^ to simple (exponential, biexponential, 
algebraic, etc.) decay laws. However, it is obvious that the NIBA gives the correct qualitative 
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picture, especially at short times. The exact dynamics is fully incoherent, yet of quite 
complicated appearance. These results again underline our earlier finding [0,11^ that the 
NIBA provides an excellent estimate for P{t) in the bulk of parameter space. The data 
shown in Fig.^ suggest that the NIBA is more accurate for K = 0.6; this indicates that the 
NIBA becomes exact for P{t) as K -^ 1/2. 

One may question how much the respective correlation functions C{t) will deviate from 
the P{t) depicted in Fig.^, since C{t) is the relevant quantity for a comparison with the 
Kondo problem. While the NIBA was shown to give an excellent approximation for the 
function P{t) even in the low-temperature Kondo region, the quality of the NIBA-prediction 
C{t) = P{t) seems to be poor in this region (see Fig.|^). We have carried out simulations 
for C{t) at temperatures slightly below Tx for 1/2 < K < 1 a.s well, and the characteristic 
behavior shown in Fig.^, namely C{t) approaching zero from below as t — i> oo, was found to 
persist. This would resolve the divergence of the static spin susceptibility |1T|J35[| predicted 
by NIBA based on P{t) because P{t) and C{t) show qualitatively different behaviors. Such 
a behavior of C{t) is also in correspondence with the exact Shiba relation |S^. 



Unfortunately, we were not able to reach the true zero-temperature limit for the equilib- 
rium correlation function C{t). Clearly, using the same Kadanoff-Baym contour employed 
in our method would require an infinitely long imaginary-time path for T = 0, which makes 
the algorithm impracticable. By invoking ergodicity arguments |jl| , however, a viable variant 
of our algorithm may facilitate such a calculation. To that purpose, one might consider a 
factorized initial state at to < 0, so that for to ~^ ~oo the system will be equilibrated at 
t = 0. In effect, one is then left with a real-time contour instead of the Kadanoff-Baym 
contour, and the initial correlations are represented by negative-time paths. This method is 
currently under study. 



21 



IV. CONCLUSIONS 

We have proposed a real-time quantum Monte Carlo simulation method for a numerically 
exact computation of the dynamical quantities of the spin-boson model. Our technique is 
based on a discretized path integral formulation and makes use of the symmetry properties 
of the dissipative influence functional, whereby one can integrate out the quasiclassical paths 
and is left with a stochastic sampling of the quantum fluctuations alone. This leads to a 
significant improvement compared to earlier methods, and allows us to study the dynamics 
at considerably longer times. The method is generally applicable to dissipative tight-binding 
systems with arbitrary spectral density. 

Results have been presented for two problems of current interest, (a) We have com- 
puted the electron transfer rate constant as a function of the electronic coupling in both the 
high-temperature and the low-temperature limit. In the high-temperature limit, the classi- 
cal Marcus result is reproduced; for low temperatures, however, one obtains an oscillatory 
behavior for large electronic couplings which does not allow for a simple rate description. 
(b) The dynamics in the antiferromagnetic Kondo region has been computed at T = for 
several values of the Kondo parameter. In accordance with the noninteracting blip approx- 
imation (NIBA), we find a fully incoherent (yet non-exponential) decay for the occupation 
probability P{t). Initial preparation effects, however, lead to important deviations from the 
NIBA-prediction C(t) = P(t) for the equilibrium correlation function. It is exactly this 
NIBA-equality which led to unphysical results for the equivalent Kondo problem. 

ACKNOWLEDGMENTS 

This work was partly supported by the Camille and Henry Dreyfus Foundation under the 
New Faculty Awards Program, the National Science Foundation (CIIE-9216221), and the 
Young Investigator Awards Program (CIIE-9257094). Computational resources furnished by 
the IBM Corporation are acknowledged. Finally, we wish to thank Uli Weiss for stimulating 



22 



discussions. 



23 



REFERENCES 

[1] U. Weiss, Quantum Dissipative Systems, Series in Modern Condensed Matter Physics, 
Vol.2 (World Scientific, Singapore, 1993), and references therein. 

[2] A.O. Caldeira and A.J. Leggett, Phys. Rev. Lett. 46, 211 (1981); Ann. Phys. (N.Y.) 
149, 374 (1983); ibid. 153, 445(E) (1983). 

[3] A.J. Leggett, S. Chakravarty, A.T. Dorsey, M.P.A. Fisher, A. Garg, and W. Zwerger, 
Rev. Mod. Phys. 59, 1 (1987). 

[4] U. Weiss, H. Grabert, and S. Linkwitz, J. Low Temp. Phys. 68, 213 (1987). 

[5] H. Wipf, D. Steinbinder, K. Neumaier, P. Gutsmiedl, A. Magerl, and A.J. Dianoux, 
Europhys. Lett. 4, 1379 (1987). 

[6] J. Kondo, Prog. Theor. Phys. 32, 37 (1964); A.M. Tsvehck and P.B. Wiegmann, Adv. 
Phys. 32, 453 (1983). 

[7] For a review, see R.A. Marcus and N. Sutin, Biochim. Biophys. Acta 811, 265 (1985), 
and references therein. 

[8] J. Ulstrup, Charge Transfer Processes in Condensed Media (Springer, New York, 1979). 

[9] D. Chandler, in Liquids, Freezing, and the Class Transition, Les Houches 51, Part 
1, edited by D. Levesque, J. P. Hansen, and J. Zinn- Justin, (Elsevier Science, North 
Holland, 1991). 

[10] We do not explicitly consider the problem of a particle in an extended double- well 
potential here. At low enough temperatures (and under certain restrictions regarding the 
system parameters) the extended problem reduces to the two-state situation discussed 
here, see, e.g., A.T. Dorsey, M.P.A. Fisher, and M.S. Wartak, Phys. Rev. A 33, 1117 
(1986). 

[11] R.P. Feynman and F.L. Vernon, Ann. Phys. (N.Y.) 24, 18 (1963). 

24 



[12] D. Chandler and P.G. Wolynes, J. Chem. Phys. 74, 4078 (1981). 

[13] C.H. Mak and D. Chandler, Phys. Rev. A 41, 5709 (1990); ibid. 44, 2352 (1991). 

[14] C.H. Mak, Phys. Rev. Lett. 68, 899 (1992). 

[15] R. Egger and U.Weiss, Z. Phys. B 89, 97 (1992). 

[16] C.H. Mak and J.N. Gehlen, Chem. Phys. Lett. 206, 130 (1993); R. Egger and C.H. 
Mak, J. Chem. Phys. 99, 2541 (1993). 

[17] V.S. Filinov, Nucl. Phys. B 271, 717 (1986). 

[18] See, for example, J.D. Doll and D.L. Freeman, in Lasers, Molecules and Methods (Adv. 
Chem. Phys., vol. LXXIII), edited by J.O. Hirschfelder, R.E. Wyatt, and R.D. Coalson 

(Wiley, New York, 1989). 

[19] N. Makri and W.H. Miller, J. Chem. Phys. 89, 2170 (1988); N. Makri, Chem. Phys. 
Lett. 193, 435 (1991). 

[20] J.E. Gubernatis, M. Jarrell, R.N. Silver, and D.S. Sivia, Phys. Rev. B 44, 6011 (1991). 

[21] R.E. Cline, Jr., and P.G. Wolynes, J. Chem. Phys. 88, 4334 (1987). 

[22] B.A. Mason and K. Hess, Phys. Rev. B 39, 5051 (1989). 

[23] E.Y. Loh, J.E. Gubernatis, R.T. Scalettar, S.R. White, D.J. Scalapino, and R.L. Sugar, 
Phys. Rev. B 41, 9301 (1990). 

[24] R. Egger and C.H. Mak, J. Phys. Chem. (in press). 

[25] C.H. Mak and R. Egger, Phys. Rev. E 49, 1997 (1994). 

[26] R. Egger, C.H. Mak, and U. Weiss, J. Chem. Phys. 100, 2651 (1994). 

[27] H. Grabert, P. Schramm, and G.-L. Ingold, Phys. Rep. 168, 115 (1988). 

[28] A. Schmid, J. Low Temp. Phys. 49, 609 (1982). 



25 



[29] The last term in the influence functional ( p.l9| ) modifies only the matrix \^^\ Note that 
the path combinatorics is already contained in the definition of the V*^-'^ matrices, since 
some of these matrix elements are zero. Finally, for a general tight-binding system with 
d states, we have to deal with a matrix multiplication of q {2d + 1 x 2d + 1) matrices. 

[30] J. P. Valleau and S.G. Whittington, in Statistical Mechanics, Part A: Equilibrium Tech- 
niques, Vol. 5 of Modern Theoretical Chemistry, ed. by B.J. Berne (Plenum, New York, 
1977). 

[31] A. Garg, J.N. Onuchic, and V. Ambegaokar, J. Chem. Phys. 83, 4491 (1985). 

[32] S.S. Skourtis, A.J.R. da Silva, W. Bialek, and J.N. Onuchic, J. Phys. Chem. 96, 8034 
(1992). 

[33] H. Grabert and U. Weiss, Phys. Rev. Lett. 54, 1605 (1985). 

[34] A. Erdelyi, Higher Transcendental Functions, Vol. 3 (McGraw-Hill, New York, 1955). 

[35] M. Sassetti and U. Weiss, Phys. Rev. A 41, 5383 (1990). 

[36] H. Shiba, Prog. Theor. Phys. 54, 967 (1975). 



26 



TABLES 
TABLE L Performance of different dynamical simulation techniques for spin-boson systems: 

(a) a brute-force evaluation of the path integrals without any filtering, (b) the SPMC calcula- 
tion [^jl^, (c) the original discrete filtering technique [14,16|, (d) the current optimized filtering 
method. The data shown below are for a calculation of P{t) for K = 1/2, ujc/ A = 6, both for 
T = and for a high temperature, A/3 = 0.5. For times t > tmax! the dynamical sign problem 
becomes uncontrollable (with statistical errors > 20%). 
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FIG. 1. Discretization of the KadanofF-Baym contour 7 in the complex-time plane. 

FIG. 2. Simulation results for K = 2, k^T/fiuOc = 4 and two values of the electronic coupling. 
Note the change in timescale. 

FIG. 3. Electron transfer rate constants as a function of the electronic coupling; squares 
denote decay constants of P{t) for X = 2, k-QT/fiijOc = 4, i.e., total rates. The dashed line is the 
nonadiabatic golden rule prediction, and vertical bars are error estimates. 

FIG. 4. Simulation results for K = 2, k-QT/fiijOc = 0.4 and two values of the electronic coupling. 
Note the change in timescale. 

FIG. 5. Symmetrized correlation function C{t) for K = 1/2, cJc/A = 6 and k^T/fil^ = 0.025. 
The triangles denote the QMC data for discretization numbers g = 60 and r = 130. The solid 
curve is the exact solution from Ref. [pH|, and vertical bars are error estimates. 



FIG. 6. Zero-temperature dynamics of P{y = Aefft) in the antiferromagnetic Kondo region. 
The solid curves are numerically exact results, dashed curves are NIB A predictions. Note the 
change in effective timescale between both plots. The discretization numbers were g = 80 and 
q = 160 for K = 0.6 and K = 0.7, respectively. 
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